function energy=Smootheddual_energy(p1,p2,f,s)
[ny,nx,n]=size(p1);
ny=ny-1;
Divp=zeros(ny,nx,n);

 Divp = p1(2:ny+1,:,:) - p1(1:ny,:,:) + p2(:,2:nx+1,:) - p2(:,1:nx,:);
 w = sum(exp( (-f - Divp)/s ),3);
 energy=s*sum(sum(log(w))); 